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Abstract 

A new generation of ultra-low-background scintillator-based detectors aims to study 
solar neutrinos and search for dark matter and new physics beyond the Standard 
Model. These optical, non-imaging detectors generally contain a "fiducial volume" 
from which data are accepted, and an "active buffer region" where there are higher 
levels of radioactive contaminants. Events are observed in real time. To distinguish 
between events occurring in the two regions, it is imperative that event position 
reconstruction be well-understood. The object of this paper is the study of the 
reconstruction, in time and space, of scintillation events in detectors of large di- 
mensions. A general, likelihood-based method of position reconstruction for this 
class of detectors is presented. The potential spatial resolution of the method is 
then evaluated. It is shown that for a spherical detector with a large number N of 
photosensitive elements that detect photons, the expected spatial resolution at the 
center of the detector is given by 5a as (ccr/n)y/3/N, where a is the width of the 
scintillator time response function and n is the index of refraction in the medium. 
However, if light in the detector has a scattering mean free path much less than the 
detector radius R, the resolution instead becomes (R/2)y / 3/N. Finally, a formalism 
is introduced to deal with the common case in which only the arrival time of the 
first photon to arrive at each photosensitive element can be measured. 
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1 Introduction 



Optical, non imaging detectors are widely used for the detection of weakly 
interacting particles. At present the main focus of observation is on neutri- 
nos and antineutrinos from various sources, but there are also plans to con- 
struct large optical detectors to search for as yet undiscovered particles such 
as WIMPs. The detection mechanism is based on the collection of visible 
or ultraviolet photons. These are emitted as Cerenkov radiation (e.g., as in 
Kamiokande [1] and SNO [2]) or as scintillation photons. We will focus our 
attention in this paper on scintillator-based, unsegmented detectors. 



1.1 A Brief History of Scintillation Detectors 

The history of scintillator-based detectors is heavily intertwined with that of 
neutrino physics. The first neutrino detector ever built, that of Cowan and 
Reines in 1953, was a 10.7ft 3 cylinder filled with a cadmium-doped organic 
scintillator and wavelength shifter, which detected reactor-generated z/ e 's by 
observing the coincidence of e + annihilation and neutron capture following 
the inverse beta decay reaction v e (p, n)e + [3,4,5]. However, the first large-scale 
unsegmented liquid scintillator detector was not built until about 1980. The 
100 ton neutrino detector at Artemovsk, Ukraine, a cylindrical 5.6 m x 5.6m 
tank filled with a saturated hydrocarbon scintillator and fluor, was a direct 
descendant of Reines and Cowan's original design. Indeed, it was designed 
to detect antineutrinos using the same reactions [6]. It was buried in a salt 
mine, 600 meters water equivalent (m.w.e.) underground, and was later used 
to study the interactions of cosmic ray muons with scintillator [7]. 

The 1995 Counting Test Facility (CTF) prototype of the Borexino exper- 
iment further developed the architecture of scintillator-based detectors [8]. 
This 4 ton detector was intended primarily as a test bed for technologies of 
the full-scale Borexino detector, not as a neutrino detector in its own right. 
Nevertheless, it set a record for the lowest detector background achieved at 
the time, of 0.03 counts/ (kg keVyr), in the window 250 keV to 2.5MeV [9]. It 
has as a result produced new upper bounds on various exotic processes [10]. 
Unlike previous scintillation detectors, it is spherical in design, in order to 
keep as much scintillator away from the surface as possible. Liquid scintilla- 
tor (both pseudocumene and phenylxylylethane, at different times, again with 
added fluors) is contained in a thin spherical nylon balloon, surrounded by 100 
inward-facing photomultiplier tubes. This setup is contained in 1000 tons of 
ultrapure water in a cylindrical tank. The entire detector is 3400 m.w.e. under- 
ground in the Gran Sasso National Laboratory, in central Italy. The CTF first 
established the feasibility of a scintillator-based solar neutrino detector with 
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a detailed study of the radioactive contaminants internal to the scintillator. It 
was also the first scintillation detector to introduce an inactive buffer (water) 
between the active volume of scintillator and the photomultiplier tubes. As 
well, it has the capability of position reconstruction for point-like events. 

The CHOOZ detector [11], built to study oscillations in reactor antineutri- 
nos from a nuclear power plant by the same name in northern France, took 
data in 1997-98. Its layered design incorporated key features of the CTF and 
Borexino designs, as well as those of other neutrino detectors such as the 
Cerenkov detector SNO [2], and the hybrid Cerenkov/scintillation light detec- 
tor LSND [12]. (The design of the larger hybrid detector MiniBooNE, built in 
1999 in order to confirm or refute results from LSND by observing 0.5-1 GeV 
muon neutrinos produced at the FNAL accelerator, was based upon the same 
principles [13].) The interior of the CHOOZ detector featured a central 5 ton 
Gd-doped target mass inside a clear roughly egg-shaped Plexiglas container, 
surrounded by an undoped 17 ton inactive buffer region contained in an ob- 
long "geode," and an outer undoped 90 ton volume with its own set of PMTs, 
used for vetoing muons from cosmic rays. The detector was placed at a depth 
of 300 m.w.e. 

The current generation of unsegmented detectors based on organic liquid scin- 
tillators - KamLAND [14], taking data since 2002, and Borexino [15], soon to 
begin operations - retain this sort of layered design, both using the spheri- 
cal shape of the CTF. Unlike the detectors described already, Borexino will 
observe scintillation light due directly to neutrino scattering from electrons, 
and can therefore potentially detect neutrinos with much lower energies (the 
threshold v energy for the inverse (3 decay is 1.8 Me V). KamLAND has ob- 
served disappearance of z/ e from reactors using the inverse f3 decay signa- 
ture [16], but it is also intended to observe solar neutrinos directly via u-e 
scattering in the future. The current KamLAND background in the region be- 
low 2 MeV must be drastically reduced for that goal to be achieved [17]. These 
detectors are situated much deeper underground (Borexino: 3400 m.w.e; Kam- 
LAND: 2700 m.w.e.), for further reduction of the residual muon flux and the 
production of short-lived cosmogenic isotopes. 

Two new experiments with targets of liquified noble gas, also aiming at low 
energy solar neutrino detection via the detection of scintillation light, are 
currently under development: CLEAN [18] and XMASS [19]. In the case of 
CLEAN, wavelength-shifter coated windows are offset from the PMTs by a 
5-10 cm gap which is a thin inactive buffer region. For these detectors, reliable 
determination of the positions of events is even more important, due to the 
need of rejecting the higher background rate coming from scintillation events 
produced in proximity of the PMTs and container vessel. An additional com- 
plication arises because the mean scattering length of scintillation photons 
(produced in the ultraviolet range of the spectrum for noble gases) is much 
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less than the radius of the detector; scintillation photons propagate from the 
event origin to the PMTs in a diffusive mode. Therefore the times of arrival 
of detected photons provide less information than in detectors using organic 
scintillator; these noble gas detectors will rely heavily upon the spatial pattern 
of PMT hits to reconstruct the positions of events. 



1.2 The Necessity of Spatial Event Reconstruction 

Due to the extremely low interaction rates of neutrinos and their antiparti- 
cles (to say nothing of WIMPs and so forth), it is necessary for a detector to 
contain a large mass of scintillator with very low levels of internal radioactive 
contamination [15]. Ultra-pure materials are also used to screen radioactivity 
from materials surrounding the detector [15,20]. Unfortunately, the photosen- 
sitive elements used to detect scintillation light are notorious for being among 
the main sources of radioactivity in an ultra-low-background detector. 

It is therefore desirable to insert, between the photosensitive elements and 
the scintillator, one or more layers of buffer material to suppress radioactive 
background. Often the buffers are inactive, i.e., not scintillating. An inactive 
buffer offers the advantage of minimizing the total trigger rate caused by 
the abundant radioactive decays generally produced within the photosensitive 
elements [15]. On the other hand, if the compositions of the scintillator and 
inactive buffer are different, a scintillator containment system is required to 
physically separate them [15]. The containment system, being in direct contact 
with the scintillator, must satisfy extremely stringent requirements in terms 
of intrinsic radiopurity. 

For additional background prevention, the outer region of the scintillator vol- 
ume can be used as an active buffer. This allows any residual radioactivity 
coming from the containment system, or passing through it, to be monitored 
and suppressed. A "fiducial volume" is commonly defined as a region at the 
center of the active volume of the detector in which radioactive background is 
expected to be at a minimum. The discrimination between events belonging 
to the fiducial and to the non-fiducial regions is performed by means of soft- 
ware implementation (reconstruction code) of an algorithm (reconstruction 
algorithm), which assigns to each single event a reconstructed position, either 
inside or outside the fiducial volume. The algorithm also provides a means of 
comparing the position of different events and is an important tool for the 
identification of several background sources. The designs of some planned de- 
tectors incorporate only a thin inactive buffer region or none at all, and in 
these cases, correct assignment of an event as belonging to the fiducial vol- 
ume or the buffer region is even more important. The resolution of detector 
reconstruction codes are generally studied with Monte Carlo methods. Event 
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simulations allow close reproductions of the performance of these codes on 
real events. Typically, however, the reconstruction codes are fine tuned by 
calibrating the detector with the use of localized sources of radioactivity or 
light. 

What seems lacking from the available literature is a comprehensive discussion 
of how the resolutions of detector reconstruction codes are related to some 
basic properties of the detector: the linear dimension, the time dispersion of 
the photon emission, the scintillator index of refraction, possible processes of 
absorption and re-emission and of scattering of the scintillator light, etc. In 
this paper we present an analytic study of the resolution for reconstruction in 
time and space of scintillation events. The study is restricted, for simplicity, to 
the case of events at the center of the detector, simple enough to be treated, 
within certain approximation, analytically. Calibrations of experiments [9] and 
full Monte Carlo studies of the performance of proposed experiments [21] show 
anyhow that the resolution of the reconstruction codes depends only in a mild 
way upon the location of the scintillation event. 

This study also assumes that the optical properties of the media are uniform 
throughout the detector, and that the indices of refraction of all materials 
between the active scintillator and the photodetectors are approximately the 
same. 



2 Likelihood Function Derivation 

The likelihood function is a standard statistical tool used for finding param- 
eters of a physical model. Suppose that a set of N observations is composed 
of the independent values {tj} and dependent values {s^} (i = 1,...,JV). For 
instance, {tj} could be a list of times at which a radioactive sample is ob- 
served, and {sj} a list of observed activities at each time. We wish to model 
the data using some function f(s) with n free parameters a. In the example, 
the function would be a decaying exponential, and the parameters would be 
the initial activity and the half-life. By definition, the likelihood function over 
the parameters is a probability distribution of obtaining the observed data 
given a specific set of parameters: 

£(a ; {(U, Si)}) = P({(t u s^} are observed | a = a ). (1) 

The difficult task is to calculate this probability based on the assumption that 
the data are correctly described by the model function f(s). Once this has 
been done, in order to calculate the most probable value of the parameters of 
the model, one simply finds the maximum of the likelihood function (or, as is 
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usually computationally easier, the minimum of — log C) in the n-dimensional 
space defined by the free parameters a. 

In the case of a scintillator-based detector, the parameters of interest are the 
position and time of an event in the detector, a = (x ,t ). The observed data 
are the positions {x{\ of the photosensitive elements, usually PMTs (inde- 
pendent values), and the times {ti} at which each element is hit by a photon 
(dependent values); i ranges from 1 to N, with iV being the number of de- 
tected photons. For now we assume that at most one photon is detected by 
each PMT, so all the distinct, and iV is also the number of PMTs that 

detect a photon. For conciseness, define the following possible events: 

• A : detector event occurs at (x ,t ) 

• B : detector hit pattern is {(x i: ti)}. 

Then, Equation (1) becomes 

£(xo,* ;{(^)}) = P(B|A). (2) 

2. 1 Factoring the Detector Likelihood Function 

Let us assume that the times at which photons are emitted by the scintillator 
are uncorrelated. Then the likelihood function will have one independent factor 
for the piece of data provided by each PMT 1 . Let the total number of working 
PMTs be T, so that N PMTs (labeled 1, . . . , N) have detected a photon, and 
T — N PMTs (labeled N + 1,...,T) have not. If we further define 

• Cj : PMT i is hit 

• Dj : PMT % detects a photon 

• Ej : PMT % detects a photon at time U, 

then 

N 

p(bia) = n p ( E *i A > c - d *) p ( d *i a > c *) p ( c *i A ) 

1=1 

x n [ p Hm Cj) P(Cj\A) + P(--C i |A)] (3) 

j=N+l 

1 Strictly speaking, this is not exactly true; specifying that N PMTs detected pho- 
tons causes the PMT hit data to be correlated. For a reasonably large number of 
hit PMTs, though, the difference should be negligible. It would be interesting to 
compare results derived from the often-used Poisson and multinomial probabilistic 
models to the model put forth here 
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(where -i is the logical negation symbol). Of course, P(Dj|A, Cj) is just the 
quantum efficiency of PMT i, which is, to a first approximation, independent 
of the original event position. 



Now define a "per-PMT" likelihood function £j. 



£i{X0, to] Xi, ti) 



= < 



5iP(Ej|A, Q, Dj) P(Q|A), i<N 
^ (1 - Qi ) P(Q|A) + P(-Q|A), N <i<T 



(4) 



The total likelihood function is then the product of all per-PMT likelihood 
functions. Notice that the per-PMT likelihood function of a supposedly dead 
PMT (qi = 0) that does not detect a photon reduces to 1, so does not influence 
the total likelihood function, just as expected. 



2.2 Scintillator Dispersion Time at the Emission Point 

The first factor in the expression for the likelihood function of a PMT that 
detects a photon is based solely on timing information of a photon emitted 
by the scintillator. Scintillation photons are emitted as a consequence of the 
ionization of the scintillator due to interacting particles or radioactive decays. 
The typical dispersion in the time of emission of organic liquid scintillators is 
on the order of a few nanoseconds, with a slower component that can reach 
hundreds of nanoseconds. The emission of photons is uniform over the solid 
angle. In this discussion we assume that the time of emission of each photon, 
relative to the time of the event causing scintillation, is an independent random 
variable r e . 

Suppose the distribution of the random variable r e is given by some scintillator 
response function p(r e ). Referring to the left half of Figure 1, one sees that at 
a specific time t, this function may also be regarded as an outgoing spherical 
photon probability wave, integrated over the solid angle 4-7T. In fact, the most 
important factor in Equation (3), the probability P(Ej|A, Q, Dj), is equal to 
it. Let rj be the time of flight from the origin xq of the photon to the position 
Xi of the i th PMT. Then, with n being the scintillator index of refraction, we 
have: 



i _ \Xi — Xq\ n 
1 c 

U = T e + Tj + t . 



(5) 
(6) 
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Fig. 1. Geometry of the likelihood function derivation. The concentric dotted lines, 
and the graph on the left, represent the probability function (an expanding spherical 
wave) of the emission time of a scintillation photon. The rectangle labeled dAi 
represents a PMT of infinitesimal size with normal vector n, subtending a solid 
angle dfij as seen from the position of the detector event. The PMT is tilted away 
from the direction of the event by an angle tpi ■ Note that we have not yet made any 
assumptions about the geometry of the detector. 

As a result, 

Ci(xo,t ;Xi,ti) oc p(ti-t -Tf). (7) 

Of course, factors other than the dispersion time of the scintillator may also 
affect the probability distribution function of the recorded arrival times of 
photons at PMTs. The most important other effects are usually the effects of 
scattering in the scintillator and the finite time resolution of the PMTs them- 
selves. The latter may in general be incorporated into the distribution p(r e ) 
by convolution with the scintillator dispersion function. The former requires a 
bit more care because scattering effects depend in general upon the light path 
length from the event to the PMT; an exact treatment is beyond the scope of 
this paper. 

2.3 Photon Attenuation 

As photons travel away from their origin, they are attenuated by the familiar 
inverse square law. This implies a formula for the probability P(Q|A) that a 
given PMT is hit by a scintillation photon. Suppose a PMT of infinitesimal 
area, at a distance Sj = \xi — xq\ from the event, subtends a solid angle dQi as 
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seen from the event location. Assuming a perfect collection efficiency, it will 
collect only a fraction df2j/47r of all photons emitted. So if T photons were 
emitted, its probability of being struck by at least one of them is 

. . . / dfL\ r dfL , . 

P(C,|A) = 1- i-^ -r^. (8) 



If the i th PMT has an area dAi and is tilted away from the line of sight by an 
angle ipi, as shown on the right half of Figure 1, then dfli = cosipi dAi/sf, so 
the resulting factor in the likelihood function is given by 

Ci{xo,to;xi,ti) oc T-— = T— — T dA^. (9) 
An Ans i 



As mentioned already, all constant factors in a likelihood function may be 
discarded with no effect on the location in parameter space of its maximum. 
(To first order, this includes the quantum efficiency qi of each PMT.) The per- 
PMT likelihood function for a PMT detecting a photon may thus be redefined 

as 

Ci(x , t ; £i, U) = p(ti - t - r}) C ° S f l . (10) 

Si 



Its logarithm is 

\ogCi = log p(ti - t - r}) + logcos^i - 21ogSi. (11) 



24 The PMTs Not Triggered 



For completeness, we now consider the case of a PMT that does not detect a 
photon produced by an event in the detector. Its per-PMT likelihood function, 
from Equation (4), is given by 



£i(xb, to) d 3 x dt 



;i-^)P(Q|A) + P(-Q|A) 

r" 



(1-%) 



1-1- 



da 



4:71 



+ 1- 



An 



= 1 - q t + q t \ \ 

da 



dfii 

An 



i - q l T- 



An 



(12) 
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« A = (aTo,fo) 



Fig. 2. Geometry of a spherical detector. 



The logarithm of this per-PMT likelihood function is — qiT&Vti/ An . This 
term, containing an infinitesimal, is negligible in size compared to the terms 
of Equation (11) coming from per-PMT likelihood functions for PMTs that 
have detected a photon. If PMTs are in fact very small compared to any other 
relevant dimensions of the detector, it may therefore be ignored. 



2.5 Specialization to a Spherical Detector 

As written, Equation (10) is applicable to any detector with pointlike PMTs 
forming the vertices of a convex polyhedron (so that light from an event at any 
point inside the detector may reach any one of the PMTs). Let us specialize 
to a spherical detector of radius R centered at the origin, having a uniform 
distribution of inward-facing PMTs over the surface. As above, we call the 
distance from an event to the i th PMT s« = \xi — xq\. Let the distance from 
the center of the detector to the event be a = \x \, so we have the geometry 
of Figure 2. 

By dropping a perpendicular from segment OB to point A (shown as line 
segment AC), one readily sees that Sicosipi = R — acos6{, with 9j being 
the angle between the event and i th PMT seen from the origin. Hence the 
likelihood function becomes 



C(x ,t ; {(x^U)}) = J[p [U - t, 



Sjn\ R — acosOi 



c ) S? 3 



(13) 



i=l 



where Sj is given by the Law of Cosines, 



s] = R 2 + a 2 - 2aR cos 6 \. 



(14) 
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3 Properties of the Likelihood Function at the Origin 



It may be of interest to examine properties of the likelihood function in the 
particular case of a hypothetical event occurring at the center of a spherical 
detector. This allows the general nature of the problem of reconstruction to 
be understood analytically. For simplicity, let's assume that the distribution 
of the time emission of the photons is a Gaussian curve with width equal to 
the characteristic dispersion time of the scintillator: 

P ^ = n\ — 9 ; lo SP( r e) = const ( 15 ) 

The same equation can also be used for the case when the original scintillation 
light is absorbed and then re-emitted by scintillation fluors in the immediate 
proximity of the energy deposition point [21]. In this case, the dispersion 
characteristic of the scintillator is effectively broadened by the absorption and 
re-emission process. 



3. 1 Taylor Expansion of the Likelihood Function 



For a point in the detector at a distance a from the center, in the direction of 
a particular unit vector u, the log likelihood function is 

1 Sin\ 2 R — acosOi 
log £(au, t ) = const - —^2^ (U - to I + 2^ log 3 (16) 

2° i=l ^ C ' i=l S i 



where Sj and 9i for each PMT shown in figure 2. We assume that 

the number of hit PMTs N is sufficiently large that we can, with little error, 
replace this expression by spatial and temporal averages over the expected an- 
gular and time distributions of the PMT hits. That is (discarding the constant 
term), 

, „, - v N If sn\ 2 \ /, R-acos9\ 
\ogC(au,t ) » -— s n t -t -—j j + Nllog -g j, (17) 

where t, s, 6 are now continuous random variables with the expected distri- 
butions. We now calculate these averages for a point-like event located in the 

— * 

center x = of the detector, occurring at time t — 0. 

First consider the time average. The time of flight of photons from the center 
to each PMT (assuming minimal scattering) is Rn/c, where n is the index 
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of refraction and c is the velocity of light in vacuum. This means that the 
distribution curve of t is p(t — Rn/c). From the properties of a Gaussian 
distribution, the time averages of time-dependent quantities are 



, . Rn , , 

(t) = — (18) 

t 2 ) = (t) 2 + a 2 = ^ + a 2 . (19) 



Likewise, since all PMTs are equidistant from an event at the center of a 
spherical detector, the distribution of PMT hits should be uniform over the 
solid angle. Hence the spatial averages over quantities dependent upon the 
event-to-PMT angle 9 can be found using Equation (14) and taking the surface 
integral over the sphere of PMTs: 



(s) = ^- /d0d(cos0) VR 2 + a 2 -2aRcos0 = R+^ (20) 

4:71 J 3R 

s 2 ^ = — J d(f) d (cos 6) (R 2 + a 2 - 2aR cos 6) = R 2 + a 2 (21) 



Finally, we observe that for a point-like event in the center of a uniform sphere 
of PMTs, there is no correlation between the expected spatial distribution of 
s and temporal distribution of t; that is, (st) = (s) (t). This and the above 
equations allow us to evaluate 



'/ sn\ 2 \ I 9 9 s 2 n 2 n n sn „ sn\ 

U -to-—) ) = (t 2 + t 2 + — - 2tt - 2t— + 2t — ) 

R 2 n 2 9 ,„ 9 9x n 2 n Rn 
= + t 2 Q + (R 2 + a 2 — - 2 — 1 

cr c z c 

Rn 2 , „ a 2 . ^ . _ a 2 .n 
= const + t 2 Q + —a 2 + ^« 2 ^o (22) 



where the constant term contains whatever does not depend explicitly on t 
and a. 

The quantity averaged over in the last term of Equation (17), again substi- 
tuting in Equation (14), becomes 
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R — a cos 9i 
log 3 = log 



R — a cos Qa 



(r 2 + a 2 - 2aRcos9, 



3/2 



2a 2 

= -21ogR + ^cos9 i + ^- s (5 cos 2 ^-3) + ... (23) 

with the last equality above being the expansion into a Taylor series in a/R. 

By once again averaging the expected distributions in s and 9 over the solid 
angle, the result, obtained to second order in a/R, is determined to be 



log 



R — a cos 9 ' 



const — 



2a^ 
3R 2 



(24) 



The complete likelihood function for an event at the center of a spherical 
detector, to second order in a/R, is thus 



log C(au, t ) w const — iV 



1 



2n 



■ (25) 



3.2 Likelihood Function Maximum and Resolutions 



Solving for the maximum of the likelihood function and requiring \a\ < R 
gives the expected solutions: 







to 
a - 



-- 




(26) 



We next ask about the expected resolution of the detector. Notice that the in- 
formation matrix is diagonal because the off-diagonal terms, — d 2 (log C) / dadt®, 
are zero when a = t = 0. The theoretical resolutions of the detector in space 
and time are therefore given by reciprocals of the second derivatives of the 
likelihood function: 




d 2 logC 

d 2 iog/^ 



-1/2 



-1/2 



a 



N 

Nn 2 
3c a 



an r 1/2 



(27) 



When the detector dimensions are much larger than the scintillator dispersion 

/ Q 

time, R ^> ca/n, we can approximate 8a ~ \ jsflT- should be noted that 
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this does not take into account scattering effects, which become increasingly 
important with larger detectors.) 

Because of the spherical symmetry of the problem, 8a can be used as a stand- 
in for any of the three Cartesian spatial resolutions Sxo, Syo, 5zo- One may, 
for instance, make the substitution a 2 = Xq + y\ + z 2 . in Equation (25) and 
obtain the same results for the resolution in each Cartesian coordinate. 



3.3 Pattern Matching 



In case of use of a liquified noble gas as scintillator, as in the new generation 
of solar neutrino detectors [18,19], Rayleigh scattering of the ultraviolet scin- 
tillation photons plays an important role. The photons are scattered intensely 
by the medium, such that they effectively diffuse out of the medium with a 
very long dispersion time; then R ^> ccr/n is no longer valid. In this case, the 
information carried by the time of flight method about the original position 
of the events becomes less reliable. However, it is still possible to reconstruct 
the original position of the event by taking into account that the density of 
hits on the PMTs decreases with the inverse of the squared distance from the 
point where the energy is deposited [22]. 

Suppose that we have no timing information, so our only information about an 
event is the pattern of hit PMTs. In this case, the likelihood function simply 
determines the position of the event. It does not depend on time and cannot 
be used to reconstruct the time itself. We may set the function p(r e ) to be 
constant and ignore it: 

„ / ^ v A. R — a cos 9i . . 

logC(au) = const + 2^ log 3 . (28) 

i=l S i 



By the same methods as above, we obtain 

log C(au) PS const —^a (29) 

3R 



for the second-order Taylor expansion in a/R of the likelihood function for an 
event at the detector center. In this case we find 

d 

— \ogC = a = 0, (30) 

(J CL 
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Table 1 

Comparison of the predicted resolutions of three liquid scintillator detectors with 
the values determined experimentally or by Monte Carlo (MC) methods. See the 
text for meanings of the columns and comments on values in italics. 

and for the resolution, 

r / dHogC Y 1 ' 2 [3R 

,a= — aH = hr (31) 



Recall Equations (25) and (27) in the case where timing information is avail- 
able: 



, 2N 2 N ( 2 n 2 2 2n , 

log £(au, t ) « const - -^a - 1 1 + + —a t 

( Nn 2 AN 

oa 




\3c z a z 3R Z J 

We see that use of timing information improves spatial resolution significantly 
when the scintillator dispersion time is much less than the travel time for 
light to cross the detector. In a liquid noble gas detector, the scintillator time 
dispersion is very broad due to the amount of internal Rayleigh scattering 
of scintillation light. Nevertheless, use of even the small amount of timing 
information available has been shown to improve the spatial resolution by a 
large fraction [23]. 



3.4 Comparison to Observed Resolutions 



Experimentally, the position resolution of a detector can be determined in 
several ways. The simplest and most common is the use of a calibration source. 
In cases when the detector has not yet been built, Monte Carlo methods are 
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of course the only method that can be used. The detector resolutions obtained 
from experimental results for CTF, and Monte Carlo tests of Borexino and 
a hypothetical liquid neon dark matter detector [23], are shown in the last 
column of Table 1. For comparison, the physical attributes of the detectors 
and the predicted resolutions 5a from Equation (27) are shown in the other 
columns of the table. As above, R is the detector radius, T the total number of 
PMTs, n the scintillator index of refraction, and a the scintillator dispersion 
time. The average number of photoelectrons detected in each event from the 
source used is denoted by e. 

iV is determined in most cases as follows. In detectors using a time-of-flight 
position reconstruction method, each PMT can measure the arrival time only 
of the first photon to strike it. This difficulty will be discussed more thoroughly 
in Section 4. The immediate consequence is that iV is a measure of the number 
of hit PMTs rather than the total number of detected photoelectrons. Basic 
probability tells us that given an event in which e photoelectrons are detected, 
the expected number of hit PMTs is 



Note, however, that for the spatial hit pattern, every photoelectron contributes 
to our knowledge, even for multiple hits on a single PMT. This implies that the 
term 4N/3R 2 in the expression for 5a in Equation (27) should in fact include 
e, not N. In calculating the predicted values of 5a in Table 1, we therefore use 
the modified expression 



Some comments on idiosyncracies of the individual detectors are in order. The 
value of n of 1.8 tabulated for the CTF is an "effective index of refraction." In 
fact, the CTF volume is partly water (n = 1.33) and partly organic scintillator 
(n = 1.5); this "effective index" is an attempt to account for refraction at the 
interface between the two fluids. Refraction causes light to travel a greater 
distance from event to PMT than it would through a single medium, so the 
"effective n" is higher than that of either pure fluid. Additionally, note that 
the observed value of 5a for the CTF takes into account only the spread in x 
and y coordinates; the CTF source had the shape of a cylinder, extended in z. 

In the hypothetical liquid neon detector described in reference [23] , events have 
a prompt component (relative intensity 2.0) and a delayed component (rela- 
tive intensity 1.0) of scintillation light. For the Monte Carlo simulation taking 
into account only the spatial pattern of PMT hits ( "spatial data only" row of 




(32) 




(33) 
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Table 1), both components contribute useful data. In that case the photoelec- 
tron yield is 2428 pe/MeV, 1.5 times the prompt light yield of 1619pe/MeV 
(1079 1.7 photons/Me V x 20% quantum efficiency x 75% geometric coverage) 
quoted in the reference. For the position reconstruction calculated from the 
spatial pattern only, we use N = e tota i = e prompt + e de i ayed in Equation (31). 

Calculation of the expected resolution in the liquid Ne detector is trickier when 
timing information is included ("timing included" row of table 1). The two 
terms contributing to 5a in Equation (33) must be evaluated with different 
values for e. The term 4e/3R 2 comes from the spatial hit pattern and so uses 
ttotai — 243, while the timing-dependent term Nn 2 /3c 2 a 2 includes only the 
prompt component of scintillation light, and thus uses e promp t = 162, with 
N = 155 derived from Equation (32). 

The source of the largest potential errors in the predictions of Table 1 is the 
value of the scintillator dispersion, a. The true scintillator dispersion function 
of a detector p(r e ) is not actually a Gaussian, so the use of Equation (27) is 
only an approximation. The value of 5.1 ns used for octf is obtained from 
the fit to CTF data described in reference [21] with the parameters shown in 
Figure 6 of that paper, sampled at 1 ns intervals and fit to a Gaussian only. 
(The same scintillator dispersion function was used in the Borexino Monte 
Carlo simulations.) Nevertheless, the predicted, observed and Monte Carlo 
values of the position resolution are in quite good agreement. For the liquid 
Ne detector, a was estimated at 10 ns, based on Figure 7 of reference [23], 
as 1/2 the difference between times with probability values equal to e -0 5 
times the value at the peak. One could plausibly estimate this value of a to 
be anywhere in the range 5.5 to 15 ns, yielding estimates of 5a from 12.6 to 
15.9 cm. This range brackets the Monte Carlo simulation nicely. 



4 Multiple PMT Occupancy and Order Statistics 

So far it has largely been assumed that the occupancy of each PMT in the 
detector is at most one. If the detector has the capability to measure the 
time at which every photon hits a given PMT, or if the detector (as with 
some of the proposed noble gas detectors) has no timing capability at all, 
then the assumption may be lifted with no effect, except that some of the Xi 
(and hence 6i and Sj) will be identical in Equation (13). For a detector with 
timing capabilities, however, it is more likely that the detector only has the 
capability to measure the arrival time of the first photon to reach each PMT. 
The probability function of the first photon to reach a PMT is not the same as 
that of a random photon reaching the same PMT; it is biased toward earlier 
times. To account for this bias, the scintillator response function p(r e ) must 
be corrected. 
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Fig. 3. A hypothetical Gaussian scintillator response function p(r e ), and its first 
order statistics for increasing values of n = 2, 3, 5, 10. Note how as n increases, the 
corrected response function narrows and shifts toward earlier times. The time axis 
is shown in units of the scintillator dispersion time a. 

4-1 Correcting for Timing Bias 

Let the probability function of the first photon to reach a PMT, out of the 
n photons reaching that PMT from an event, be represented by p n (r e ). This 
is known as the "first order statistic." Naturally, pi(r e ) = p(r e ). In general, 
the corrected scintillator response function p corr would then be some linear 
combination of the first order statistics, 



and an a priori guess would have to be made for the probability that each 
possible number of photons had hit the PMT. For simplicity, let us assume 
that the number of photons striking each PMT for an event is known (in 
Borexino, for instance, this is determined via ADC channels separate from 
the timing channels). We can then set p corr equal to the function p n {r e ). 

It remains only to calculate p n {r e ) given p{r e ) and n. Number the emission 
time of the n photons detected by a given PMT in some specific but randomly 
chosen order (for instance, in order of increasing longitude of their emission 
directions), T\, . . . , r n . Also number them in order of increasing emission time, 
si,...,s n . Then p n {r e ) is the probability function of the randomly chosen 
emission time T\ given that S\ — T\. 



■DC' 



Pcorri'Te) = ^2 Pn( T e) x P(n photons hit the PMT), 



(34) 



n=l 
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p n (r e )dr e = P(ti G [r,r + dr]|ri = si) 

= P(ti = si|ti G [t, r + dr]) x 



P(ri G [r,r + drj) 

P(n = si) 



p(T e ) dr e 
(1/n) 



P(ti = Si|t! G [r,r + dr]), 



(35) 



where the second equality is once again due to Bayes' Theorem. The proba- 
bility in the last line above is just the probability that every other photon has 
a later arrival time than the randomly selected value t\. 



P(ti = si\n G [r, r + dr]) = J] P(n > n|ri G [r, r + dr]) 



i=2 



= P(r 2 > Ti|ti G [r,r + dr]) 



n-l 



/ P(0 dr^ 



n-l 



(36) 



Hence (letting -F(r e ) = /I^ p(r e )dr e represent the cumulative distribution 
function of r e ), the first order statistic of p(r e ), if n photons are detected by a 
given PMT, is 



Pn(r e ) = np(r e ) [1 - F(r e )] 



n-l 



(37) 



Graphs of the first order statistics of a representative scintillator response 
function are shown in Figure 3 for values of n equal to 1, 2, 3, 5, and 10. (The 
specific response function shown is a Gaussian, Equation (15) offset by five 
units of a from time zero.) Note how as n increases, the time distribution of 
the first PMT hit narrows and shifts toward earlier times. 



4-2 Effects on Detector Resolution 



One may ask about the effect of this correction on the likelihood function 
and spatial resolution. Consider again the case of a Gaussian scintillator time 
response function. We have 

logp„(r e ) = const + logp(r e ) + (n - 1) log[l - F(r e )]. (38) 

Substituting in F{r e ) = (1 + erf(r e /er\/2))/2, the Taylor expansion to second 
order in r e becomes 

logp„(r e ) = const - („ - - (\ + — ) % + 0(r e 3 ). (39) 
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That is, the first photon detected at each PMT contributes to the log of the 
likelihood function in the amount of —r 2 /2a 2 , but each additional photon 
contributes only in the amount of —T 2 /na 2 (plus a term linear in r e which 
has relatively little effect on the resolution for a large detector); compare 
to Equation (15). The resolution is better than if the corrected scintillator 
response function were not used, but still poorer than if the time of arrival of 
every photon could be measured. 

Suppose that the total number of photons detected is e, by N PMTs, and in 
particular that the i th PMT sees n« photons. Denoting the emission time by 
t* = ti — t — SiTi/c, the general likelihood function is then 



log C(au, t ) = const - £ Q + - L ^— j ( r e) 



1 N ,■ R-acosOi , . 

]T(n, - 1) rl + £ log -3 >-. (40) 



Define the excess photon multiplicity as 5 = (e—N)/N. The likelihood function 
in the limit of homogeneous PMT coverage as iV — > oo, for an event at the 
detector center, becomes 



log£(cra,to) = const - ^- + - J ^(t*) ^ 

N6 [2 , ,\ Ar/r _ /, R-acos9j\ 

-vv^> + iV( * +1) \ log ^r^)- 

Running through calculations analogous to those of Section 3.1, we finally ob- 
tain the explicit function 



N I 5\ ( o n J 2 2n , 
logC{au,t ) = const - -j I - + - I 1 1 + + —a t 




-N{5 + l)— 2 a\ (41) 



In the limit ca/R — > (that is, for a very large detector compared to the 
width of the scintillator response function), it can be shown that the spatial 
resolution at the center of a detector, with N and 5 varying while holding e 
constant, is proportional to \Jtt(1 + 5)/^/ti + 25. Hence the resolution of an 
event with an average photon multiplicity of 5 = 0.5 excess photons per PMT 
is 6.7% worse than if PMTs could detect the arrival time of every photon. 
With 5 = 1 excess photon per PMT (every hit PMT seeing an average of 2 
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photons), the resolution is 10.5% worse. In the limit of large 5 (for instance 



with a high-energy event), the resolution reaches an asymptote of y tt/2 times 
(about 25.3% worse) that of an ideal detector observing an event of equal 
energy. 

Realistically, construction of an ideal detector, one that measures the time of 
arrival for every photon, would be non-trivial. One may on the other hand ask, 
given a detector capable of measuring time of arrival only for the first photon 
at each PMT, how the use of the statistically corrected scintillator dispersion 
function improves the results over the use of an uncorrected function. This 
comparison is equivalent to fixing N while (for the uncorrected dispersion 
function) setting 5 to zero. In this case, the use of the corrected dispersion 
function is an improvement by the factor s/tt/ y/n + 25 (recall that smaller 
resolutions are better). For 5 = 0.5, the reciprocal of the improvement factor 
is 1.15, and for 5 — 1, it is 1.28; for large 5, it would theoretically improve 
without bound. This analysis even leaves aside the fact that for events offset 
from the center of the detector, use of the uncorrected scintillator dispersion 
function will produce a statistically biased position estimate. 



5 Conclusions 

We analyzed the resolution of spherical, optical, non-imaging scintillation 
based detectors in reconstructing the position of point-like events, limiting 
the analytic derivation to the case of events near the center of the detector. 
We found that the fundamental length scale of the resolution given by the 
time of flight method is proportional to the product of the speed of light 
in the medium and the dispersion time at the scintillation emission, as in 



In case the dispersion of the scintillation photons arrival times grows above 
the ratio of the speed of light to the detector radius, the time of flight method 
no longer gives relevant information about the point of origin of the event. The 
position of the event can still be determined by the analysis of the density of 
hits, and in this case the fundamental resolution is set by the radius of the 

/ O D 

detector, as in 5a = \ jrrQ- 



Finally, we made some comments on the need to correct the scintillation dis- 
persion function in the common case where PMT hit timing information is 
only available for the first photon to strike each PMT. In this case, even with 
a corrected scintillation dispersion function, the spatial resolution will be up 
to 25 percent worse for high-energy events compared to a similar detector 
capable of measuring timing information for all photons. 
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